home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48hor1 / factor.doc < prev    next >
Text File  |  1995-03-31  |  7KB  |  208 lines
  1. (Comp.sys.handhelds) 
  2. Item: 1892 by jurjen at cwi.nl 
  3. Author: [Jurjen NE Bos] 
  4.   Subj: HP48: The ultimate factoring program (up to 18 digits) 
  5.   Date: Fri Feb 01 1991 16:22  
  6.  
  7. This is a program that can factor large integers (up to 63 bits) on the 
  8. HP48.  It is inspired by the postings of Klaus Kalb of last week.  In fact, I 
  9. used two of Klaus' routines (see below). 
  10. The program has the following features: 
  11. - very fast (average 4 seconds for numbers of 12 digits, 
  12.   1 minute for 18 digits.) 
  13. - uses two different algorithms for factoring (trial division and Pollard rho) 
  14.   and proves all factors prime (using Selfridge's theorem). 
  15. - combines RPL (the language in the manual) with forth (what's in the ROM) and 
  16.   machine code.  This gives speed and servicability at the same time. 
  17. - The recursive calls are nicely shown during the computation. 
  18. - A separate prime testing routine. 
  19.  
  20. There's also a bug: 
  21. - It doesn't quite work with wordsizes that are a multiple of four.  This is 
  22.   due to the carries.  If you set the wordsize to 64 bits, it will be 
  23.   modified to 63 to make it work. 
  24.  
  25. [Note: Jurgen fixed this bug; see his posting at the end.  I put his bugfix 
  26. into his program; this debugged version is what's in FACTOR on the disk. -jkh-] 
  27.  
  28. How to use: 
  29. Download the object at the end of this posting.  Use ASC-> to convert it into 
  30. a directory.  Save it under the name FACTOR.  Enter the directory.  Press 
  31. DEMO.  See what happens. 
  32. Remark: factoring numbers with large factors (more than 7 digits) can take 
  33. up to ten minutes. 
  34. Mind the bug above: if you set the wordsize to a multiple of four, factoring 
  35. might take forever.  [No, it won't; see bugfix note above.  -jkh-] 
  36.  
  37. Some interesting numbers to try: 
  38. 3^37+1 (to get this, type 63 STWS #3 #37 #1 NEG POWMOD 1 +).  The program 
  39.   thinks it found a large prime factor, but quickly finds out it isn't. 
  40. 100264053529.  This number looks even more prime. 
  41.  
  42. (Intermezzo:) 
  43. Some puzzles for experienced HP48 users: 
  44.         (For information, send Email.  Please realize these are PUZZLES, so 
  45.         they aren't trivial, and explanation is complicated.) 
  46. - Compute 'SIN(180_\^o)' .  You do not get 0.  Why? 
  47. - The time to compute '253!' is much longer than the time to compute '253.1!'. 
  48. - To compute the factorial for large negative x, use 
  49. '\pi*x/(SIN(\pi*x)*(-x)!)' 
  50.   instead of 'x!' for faster results.  (Don't forget RAD.) 
  51.  
  52. Happy Hacking, 
  53.         Jurjen Bos 
  54.  
  55. -------------------Listing of the FACTOR directory (do not download)------ 
  56. For the nosy people among you, here the listing of the routines: 
  57. [FYI; download FACTOR or FACTOR.ASC from disk.  -jkh-] 
  58.  
  59. DEMO    @ Factor a random integer and time the computation 
  60. \<< DEC TICKS RAND 1000000 * R\->B RAND 1.E12 * + RAND 1.E18 * + FACTOR 
  61.   TICKS ROT - B\->R 1.220703125E-4_s * SWAP 
  62. \>> 
  63.  
  64. PRIME?  @Prime test     [NOTE: See bugfix posting below.  -jkh-] 
  65. \<< # 0d + 
  66.   CASE DUP # 2d < 
  67.     THEN DROP 0 
  68.     END DUP # 210d BGCD # 1d == 
  69.     THEN RCWS 63 MIN STWS 0 'L' STO LPRT 'L' PURGE 
  70.     END DUP # 121d < 
  71.     THEN B\->R { 2 3 5 7 } SWAP POS 0 \=/ 
  72.     END DROP 0 
  73.   END 
  74. \>> 
  75.  
  76. FACTOR  @The factoring routine  [NOTE: See bugfix posting below.  -jkh-] 
  77. \<< RCWS 63 MIN STWS # 0d + 0 'L' STO LFAC 'L' PURGE 
  78. \>> 
  79.  
  80. MULMOD  @Modulo multiplication of binaries 
  81. (In machine code. 
  82. a b n -> a*b mod n 
  83. ) 
  84.  
  85. POWMOD  @Modulo powering of binaries 
  86. (In machine code. 
  87. a e n -> a^e mod n 
  88. ) 
  89.  
  90. BGCD    @Greatest commond divisor of binaries 
  91. (In machine code.  Written by Klaus Kalb. 
  92. ) 
  93.  
  94. TRIAL   @Trial division routine 
  95. (In machine code and forth.  Machine code part written by Klaus Kalb. 
  96. ) 
  97.  
  98. MILRAB  @Miller-Rabin test 
  99. (In forth. 
  100. n base -> 1 if n is pseudoprime  
  101. ) 
  102.  
  103. LFAC    @Internal factoring routine 
  104. \<< DUP 'L' INCR DISP "Trial division" 'L' INCR DISP # 2000d TRIAL 
  105.   IF DUP # 1d == 
  106.   THEN DROP { } 
  107.   ELSE 1 \->LIST        @List of composite factors 
  108.   END 
  109.   WHILE DUP SIZE 
  110.   REPEAT SWAP OVER 1 GET 
  111.     CASE DUP # 4012009d < 
  112.       THEN B\->R + 
  113.       END DUP L 1 - DISP "Primetest" L DISP DUP LPRT 
  114.       THEN 
  115.         IF DUP # 1000000000000d < 
  116.         THEN B\->R 
  117.         END + 
  118.       END "Pollard \Gr" L DISP 
  119.       WHILE P\Gr DUP # 1d == 
  120.       REPEAT DROP 
  121.       END 2 \->LIST ROT SWAP + SWAP 
  122.     END SWAP 2 OVER SIZE SUB 
  123.   END DROP 'L' " 
  124. " OVER DECR DISP 1 STO- 
  125. \>> 
  126.  
  127. LPRT    @Internal primtesting routine 
  128. \<< 
  129.   CASE DUP # 3d MILRAB NOT 
  130.     THEN 0 
  131.     END DUP # 25000000000d > 
  132.     THEN SRT DUP 
  133.     END DUP # 2d MILRAB NOT 
  134.     THEN 0 
  135.     END DUP # 1373653d < 
  136.     THEN 1 
  137.     END DUP # 5d MILRAB NOT 
  138.     THEN 0 
  139.     END DUP # 25326001d < 
  140.     THEN 1 
  141.     END DUP # 7d MILRAB NOT 
  142.     THEN 0 
  143.     END DUP # 3215031751d \=/ 
  144.   END SWAP DROP 
  145. \>> 
  146.  
  147. SRT     @Test primality of large numbers with Selfridge test. 
  148. \<< DUP 1 - LFAC # 3d \-> n l a 
  149.   \<< "Selfridge" L DISP 9 CF 1 1 l SIZE 
  150.     FOR k 
  151.       IF l k GET SWAP OVER \=/ 8 FC?C OR 
  152.       THEN 
  153.         CASE a n 3 PICK / n POWMOD # 1d DUP2 \=/ 
  154.           THEN SWAP 3 PICK # 0d + n POWMOD \=/ 8 + SF 
  155.           END + 'a' STO+ RAND .1 \>= 
  156.           THEN 
  157.           END n a MILRAB 
  158.           THEN 
  159.           END 9 SF 
  160.         END 
  161.       END 8 FS? 9 FS? - 
  162.     STEP DROP 9 FC?C 
  163.   \>> 
  164. \>> 
  165.  
  166. P\Gr    @Pollard rho factoring algorithm. 
  167. \<< \-> n 
  168.   \<< RAND n B\->R * R\->B DUP NEWOB 
  169.     WHILE n Code n BGCD DUP # 1d == 
  170.     @ Code is a machine code routine that does 50 rounds 
  171.     REPEAT DROP 
  172.     END ROT ROT DROP2 n 
  173.   \>> OVER / 
  174. \>> 
  175. --------------------------------- end of FACTOR directory -------------------- 
  176.  
  177. (Comp.sys.handhelds) 
  178. Item: 1893 by jurjen at cwi.nl 
  179. Author: [Jurjen NE Bos] 
  180.   Subj: HP48: factoring program, patch 
  181.   Date: Fri Feb 01 1991 16:22  
  182.  
  183. The bug mentioned in the previous posting can be succesfully solved by the 
  184. following minor patches.  This is typed by hand.  The easiest way to patch is 
  185. just using VISIT instead of trying to download it. 
  186.  
  187. Happy hacking, 
  188.         Jurjen Bos 
  189.  
  190. [NOTE: The following bugfix has already been incorporated into the version of 
  191. FACTOR that's on this disk.  It is listed here FYI only.  -jkh-] 
  192.  
  193. PRIME? 
  194. \<< # 0d + # FFFFFFFFFFFFFFFFh SR AND 
  195.   CASE DUP # 2d < 
  196.     THEN DROP 0 
  197.     END DUP # 210d BGCD # 1d == 
  198.     THEN 0 'L' STO LPRT 'L' PURGE 
  199.     END DUP # 121d < 
  200.     THEN B\->R { 2 3 5 7 } SWAP POS  0 \=/ 
  201.     END DROP 0 
  202.   END 
  203. \>> 
  204.  
  205. FACTOR 
  206. \<< # 0d + # FFFFFFFFFFFFFFFFh SR AND 0 'L' STO LFAC 'L' PURGE 
  207. \>> 
  208.